#include "StdAfx.h"

// ftp://quasar.ipa.nw.ru/incoming/EPM2008/ 

#include "EphemEPM.h"
#include "CustomException.h"

double CEphemEPM::s_dRho = 82.3005676536174207; //const Earth/Moon + 1

CEphemEPM::CSingleBody::CSingleBody()
: m_dJdt(0.0)
, m_dDelta(0.0)
, m_nCoef(0)	
, m_nDim(0)		
, m_nBl(0)
{
}

bool CEphemEPM::CSingleBody::InitTxt(LPCTSTR stzTxtFilePath)
{
	if(m_nCoef != 0)
		throw CErrorAlreadyInit();

	FILE *lfn = fopen(stzTxtFilePath, _T("rb"));
	if( lfn == NULL ) 
		throw CErrorNoFile(stzTxtFilePath);
		//throw CCustomException(_T("CEphemEPM::CSingleBody::InitTxt: can not open file %s"), stzTxtFilePath);

	int nJd = 0;
	fscanf(lfn, _T("%i"),	&nJd);
	fscanf(lfn, _T("%lE"),	&m_dJdt);
	fscanf(lfn, _T("%lE"),	&m_dDelta);
	fscanf(lfn, _T("%i"),	&m_nCoef);
	fscanf(lfn, _T("%i"),	&m_nDim);
	fscanf(lfn, _T("%i"),	&m_nBl);

	m_dJdt += nJd;

	int coef_size = m_nCoef * m_nDim * m_nBl;
	m_vecCoef.resize(coef_size);

	double *pCoef = &m_vecCoef[0];
	double *pEndBuff = pCoef + coef_size;

	for( ; pCoef < pEndBuff; ++pCoef)
		fscanf(lfn, _T("%lE"), pCoef++);
	

	fclose(lfn); 
	return true;
}

bool CEphemEPM::CSingleBody::InitBin(LPCTSTR stzBinFilePath)
{
	if(m_nCoef != 0)
		throw CErrorAlreadyInit();

	FILE *lfn = fopen(stzBinFilePath, _T("rb"));
	if(lfn == NULL) 
		throw CErrorNoFile(stzBinFilePath);
		//throw CCustomException(_T("CEphemEPM::CSingleBody::InitTxt: can not open file %s"), stzBinFilePath);

	int nJd = 0;
	fread(&nJd,		sizeof(nJd),		1, lfn);
	fread(&m_dJdt,	sizeof(m_dJdt),		1, lfn);
	fread(&m_dDelta,sizeof(m_dDelta),	1, lfn);
	fread(&m_nCoef, sizeof(m_nCoef),	1, lfn);
	fread(&m_nDim,	sizeof(m_nDim),		1, lfn);
	fread(&m_nBl,	sizeof(m_nBl),		1, lfn);
	
	m_dJdt += nJd;

	int coef_size = m_nCoef * m_nDim * m_nBl;
	m_vecCoef.resize(coef_size);
	fread(&m_vecCoef[0], sizeof(m_vecCoef[0]), coef_size, lfn);

	fclose(lfn); 
	return true;
}

double CEphemEPM::CSingleBody::GetMass() const
{
	return m_dMass;
}

void CEphemEPM::CSingleBody::SetMass(double dMass)
{
	m_dMass = dMass;
}

int CEphemEPM::CSingleBody::calc(double dJdt, double* xv, double *vv, double *dMass /*= NULL*/, bool bGetVelocity/* = true*/)
// Input
//    int dJdt          - Julian date at which interpolation is wanted
// Output
//    double* xv      - object positions
//    double* vv      - object velocities
{ 
	if(m_nCoef == 0)
		throw CErrorNotInit();

	double td;             // - Julian in [initial date, final date] -> td in [-1, 1]
	int cur_bl;            // - number of approximation interval for current Julian date     
	double tt[32], xx[32]; // - arrays of values and integrals of Chebyshev polynomials
	int i;                 // - counter

	td = (dJdt - m_dJdt) / m_dDelta;

	cur_bl = (int)td;
	td = td - cur_bl;

	td = 2.0 * td - 1.0;

	if((cur_bl < 0) || (cur_bl >= m_nBl)) 
		throw CErrorOutOfRange(dJdt);
		//throw CCustomException(_T("Exceed ranges, cur_bl = %d must be in [0, %d]"), cur_bl, m_nBl);

	double *pCoef = &m_vecCoef[cur_bl * m_nCoef * m_nDim];

	calc_cheb_pol(td, tt, xx);    

	for(i = 0; i < m_nDim; i++)
	{ 
		double * pCurrentCoef = pCoef + m_nCoef * i;

		xv[i] = calc_pos(pCurrentCoef, xx, m_dDelta);
		if(bGetVelocity)
			vv[i] = calc_vel(pCurrentCoef, tt);
	}

	if(dMass)
		*dMass = m_dMass;

	return true;
}

double CEphemEPM::CSingleBody::calc_vel(double *pCoef, double* tt)
{
	double v = 0.0;
	for(int i = (m_nCoef-2); i > 0; i--)
		v = v + pCoef[i] * tt[i];
	v = v + pCoef[0];
	return v;
}

double CEphemEPM::CSingleBody::calc_pos(double *pCoef, double* xx, double dlt)
{
	int i;
	double x = 0.0;
	for (i = (m_nCoef-2); i >= 0; i--) 
		 x = x + pCoef[i] * xx[i]; 

	x = 0.5 * dlt * x + pCoef[(m_nCoef-1)];
	return x;
}

void CEphemEPM::CSingleBody::calc_cheb_pol(double td, double *tt, double *xx)
{       
	int i, j;                // counters
	bool flag;                  // (-1)^n

	double d = 2.0 * td; 
	
	{
		double * ptt_0 = tt + 0;
		double * ptt_1 = tt + 1;
		double * ptt_2 = tt + 2;

		* ptt_0 = 1.0;
		* ptt_1 = td;

		for (i = 2; i < m_nCoef; ++i, ++ptt_0, ++ptt_1, ++ptt_2) 
			*ptt_2 = *ptt_1 * d - *ptt_0;
			//tt[i] = tt[i - 1] * d - tt[i - 2]; 
	}

	{
		double * ptt_0 = tt + 3;
		double * ptt_2 = tt + 1;

		xx[0] = td;
		xx[1] = (tt[2] + tt[0]) * 0.25;
		double * pxx = xx + 2;

		for (i = 3; i < m_nCoef; ++i, ++pxx, ++ptt_0, ++ptt_2) 
			*pxx  = 0.5 * ( *ptt_0 / i - *ptt_2 / (i-2)); 
	}

	j = 0; 
	flag = false; 
	i = 3; 

	while (i < m_nCoef-1) 
	{ 
		j = j + 1; 
		d = 0.25 / j + 0.25 / (j + 1); 
		flag = !flag; 
		if (flag) 
			d = -d; 
		xx[i] = xx[i] + d; 
		i = i + 2; 
	}
}

// =============================================== //

static bool CheckFile(LPCTSTR filename)
{
	FILE * f = fopen(filename, "r");
	fclose(f);
	return (f != 0);
}

CEphemEPM::CEphemEPM(LPCTSTR stzDataFolderPath)
{
	struct  
	{
		eBodies body;
		LPCTSTR stzBobyName;
		LPCTSTR stzBinFileName;
		LPCTSTR stzTxtFileName;
		double  dMass;
	} arrFileNames [] = {													  		
		eMERCURY,	_T("mercury"),	_T("mercury.08b"),	_T("mercury.08t"),        49125.4957e-15,
		eVENUS,		_T("venus"),	_T("venus.08b"),	_T("venus.08t"),	     724345.2333e-15,
		eEARTH,		_T("earth_m"),	_T("earth_m.08b"),	_T("earth_m.08t"),	     899701.1347e-15,
		eMARS,		_T("mars"),		_T("mars.08b"),		_T("mars.08t"),		      95495.4869e-15,
		eJUPITER,	_T("jupiter"),	_T("jupiter.08b"),	_T("jupiter.08t"),	  282534584.0787e-15,
		eSATURN,	_T("saturn"),	_T("saturn.08b"),	_T("saturn.08t"),	   84597060.7310e-15,
		eURANUS,	_T("uranus"),	_T("uranus.08b"),	_T("uranus.08t"),	   12920248.2576e-15,
		eNEPTUNE,	_T("neptune"),	_T("neptune.08b"),	_T("neptune.08t"),	   15243591.0921e-15,
		ePLUTO,		_T("pluto"),	_T("pluto.08b"),	_T("pluto.08t"),	       2166.8443e-15,
		eMOON,		_T("moon"),		_T("moon.08b"),		_T("moon.08t"),				  0.0,
		eSUN,		_T("sun"),		_T("sun.08b"),		_T("sun.08t"),		0.01720209895*0.01720209895,
	};

	
	arrFileNames[eMOON].dMass  = arrFileNames[eEARTH].dMass / s_dRho;
	arrFileNames[eEARTH].dMass = arrFileNames[eEARTH].dMass - arrFileNames[eMOON].dMass;

	TCHAR stzPath[MAX_PATH];
	int nPathLen = strlen(stzDataFolderPath);
	memcpy(stzPath, stzDataFolderPath, nPathLen * sizeof(TCHAR));
		
	m_vecBodyes.resize(sizeof(arrFileNames)/sizeof(arrFileNames[0]));
	_ASSERT(m_vecBodyes.size() == BODY_COUNT);

	for(size_t pos = 0; pos < m_vecBodyes.size(); ++pos)
	{
		if( arrFileNames[pos].body != pos )
			throw CError("UnAvailable bodies");

		bool bDone = false;
		for(int nBin = 0; nBin <= 1; ++nBin)
		{
			LPCTSTR stzFileName = nBin? arrFileNames[pos].stzBinFileName: arrFileNames[pos].stzTxtFileName;
			sprintf(stzPath + nPathLen, _T("/%s"), stzFileName);
			if(CheckFile(stzPath))
			{
				if(nBin)
					m_vecBodyes[pos].InitBin(stzPath);
				else
					m_vecBodyes[pos].InitTxt(stzPath);

				bDone = true;
				break;
			}
		}
		m_vecBodyes[pos].SetMass(arrFileNames[pos].dMass);
		if(!bDone) 
			throw CErrorNoFile(arrFileNames[pos].stzBobyName);
			//throw CCustomException(_T("can not load %s data"), arrFileNames[pos].stzBobyName);
	}
}

int CEphemEPM::calc_bar(int pl_num, double dJnd, double* xv, double* vv)
//Calculates barycentric coordinates and velocities of the following objects
// 1- Mercury, 2- Venus  , 3- Earth,  4- Mars,  5- Jupiter 
// 6- Saturn,  7- Uranus , 8- Neptune, 9 - Pluto,     
// 10- Moon, 11- Sun, 13 - Earth + Moon 
// Input
//    int pl_num      - number of object
//    int jd          - Julian date (integral part) at which interpolation is wanted
//    double dj       - Julian date (frac part)
// Output
//    double* xv      - object positions
//    double* vv      - object velocities
//    int status      - status of output data
{
	pl_num --;

	double xv_add[3];      // - additional positions
	double vv_add[3];      // - additional velocities  
	int i;

	if(pl_num == eEARTH) //Earth
	{//Earth(bar) = (Earth + Moon)(bar) - Moon(geo)/rho
		if(!m_vecBodyes[eEARTH].calc(dJnd, xv, vv)) 
			return 0;
		if(!m_vecBodyes[eMOON].calc(dJnd, xv_add, vv_add)) 
			return 0;

		for(i = 0; i < 3; i++)
		{
			xv_add[i] = xv_add[i] / s_dRho;
			vv_add[i] = vv_add[i] / s_dRho;
			xv[i] = xv[i] - xv_add[i];
			vv[i] = vv[i] - vv_add[i];       
		}
		return 1;
	}

	if(pl_num == eMOON) //Moon
	{//Moon(bar) = Earth(bar) + Moon(geo)
		if(!m_vecBodyes[eEARTH].calc(dJnd, xv, vv)) 
			return 0;
		if(!m_vecBodyes[eMOON].calc(dJnd, xv_add, vv_add)) 
			return 0;

		for(i = 0; i < 3; i++)
		{
			xv_add[i] = xv_add[i] / s_dRho;
			xv[i] = xv[i] - xv_add[i];
			vv_add[i] = vv_add[i] / s_dRho;
			vv[i] = vv[i] - vv_add[i];       
		}

		if(!m_vecBodyes[eMOON].calc(dJnd, xv_add, vv_add)) 
			return 0;
		for(i = 0; i < 3; i++)
		{
			xv[i] = xv[i] + xv_add[i];
			vv[i] = vv[i] + vv_add[i];       
		}
		return 1;
	}
	if(pl_num == 12) //Earth + Moon
		return m_vecBodyes[eEARTH].calc(dJnd, xv, vv);

	// neither Earth, nor Moon
	return m_vecBodyes[pl_num].calc(dJnd, xv, vv);
}

int CEphemEPM::calc_EPM(int pl_num, int centr_num, double dj, double* xv, double* vv)
//Calculates coordinates and velocities of the oblect 'PL_NUM'
//with respect to the point 'CENTR_NUM'
//The numbering convention for 'PL_NUM' and 'CENTR_NUM':
//                                                                               
//                1 = MERCURY           8 = NEPTUNE                              
//                2 = VENUS             9 = PLUTO                                
//                3 = EARTH            10 = MOON                                 
//                4 = MARS             11 = SUN                                  
//                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER              
//                6 = SATURN           13 = EARTH-MOON BARYCENTER                
//                7 = URANUS                               
// Input
//    int pl_num      - number of object
//    int centr_num   - number of object
//    int jd          - Julian date (integral part) at which interpolation is wanted
//    double dj       - Julian date (frac part)
// Output
//    double* xv      - object positions
//    double* vv      - object velocities
//    int status      - status of output data
{
	double xv_bar[3];      // - barycentric positions of the point 'CENTR_NUM'
	double vv_bar[3];      // -  -"- velocities  -"-

	memset(xv, 0, sizeof(xv_bar) );
	memset(vv, 0, sizeof(vv_bar) );

	if((pl_num == centr_num) || (pl_num == 12)) 
		return 1;

	if((centr_num == 3) && (pl_num == 10)) // Moon(geo)
		return m_vecBodyes[pl_num-1].calc(dj, xv, vv);

	if(centr_num == 12) //barycentric
		return calc_bar(pl_num, dj, xv, vv);

	// general case 
	if(!calc_bar(pl_num, dj, xv, vv)) 
		return 0;
	if(!calc_bar(centr_num, dj, xv_bar, vv_bar)) 
		return 0;

	for(int i = 0; i < 3; i++)
	{
		xv[i] = xv[i] - xv_bar[i];
		vv[i] = vv[i] - vv_bar[i];       
	}
	return 1;
}

bool CEphemEPM::GetCoord(double date, bool bBaryCenter, AllBodies & coords, bool bGetVelocity)
{
	memset(&coords, 0, sizeof(coords));

	for(size_t bodyNum = 0; bodyNum < BODY_COUNT; ++bodyNum)
	{
		SBody & bodyRef = coords.bodies[bodyNum];
		m_vecBodyes[bodyNum].calc(date, bodyRef.x, bodyRef.v, &bodyRef.GM, bGetVelocity);
	}

	for(int i = 0; i < 3; i++)
	{
		coords.bodies[eEARTH].x[i] -= coords.bodies[eMOON].x[i] / s_dRho;
		coords.bodies[eEARTH].v[i] -= coords.bodies[eMOON].v[i] / s_dRho;

		coords.bodies[eMOON].x[i] += coords.bodies[eEARTH].x[i];
		coords.bodies[eMOON].v[i] += coords.bodies[eEARTH].v[i];
	}

	if(!bBaryCenter)
	{
		SBody & sun = coords.bodies[eSUN];
		for(size_t bodyNum = 0; bodyNum < BODY_COUNT; ++bodyNum)
		{
			SBody & bodyRef = coords.bodies[bodyNum];
			for(int i = 0; i < 3; i++)
			{
				bodyRef.x[i] -= sun.x[i];
				bodyRef.v[i] -= sun.v[i];
			}
		}
	}
	return true;
}

